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Abstract Using the dedicated computer Janus, we follow the nonequilibrium dynamics of 
the Ising spin glass in three dimensions for eleven orders of magnitude. The use of integral 
estimators for the coherence and correlation lengths allows us to study dynamic hetero- 
geneities and the presence of a replicon mode and to obtain safe bounds on the Edwards- 
Anderson order parameter below the critical temperature. We obtain good agreement with 
experimental determinations of the temperature-dependent decay exponents for the ther- 
moremanent magnetization. This magnitude is observed to scale with the much harder to 
measure coherence length, a potentially useful result for experimentalists. The exponents 
for energy relaxation display a linear dependence on temperature and reasonable extrapola- 
tions to the critical point. We conclude examining the time growth of the coherence length, 
with a comparison of critical and activated dynamics. 
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1 Introduction 

Below their glass temperature, Spin Glasses [1] (SG) are perennially out of equilibrium. The 
understanding of their sophisticated dynamical behavior is a long standing challenge both 
to theoretical and to experimental physics. 

Aging 1 2 1 is a feature of SG dynamics that shows up even in the simplest experimental 
protocol, the direct quench. In these experiments, the SG is cooled as fast as possible to 
the working temperature below the critical one, T < T c . It is then let to equilibrate for a 
waiting time, f w , its properties to be probed at a later time, t + f w . For instance one may 
cool the SG in the presence of an external field, which is switched off at time f w . The so- 
called thermoremanent magnetization decays with time, but the larger ? w is, the slower the 
decay. In fact, it has been claimed that, if the cooling is fast enough, the thermoremanent 
magnetization depends upon t and r w only through the combination f/f w , at least for 10~ 3 < 
r/r w < 10 and ? w in the range 50 s — 10 4 s 151 . In other words, the only characteristic time 
scale is the sample's own age as a glass, t„ (this behavior is named Full Aging). Note, 
however, that there is some controversy regarding the natural time variable which could 
rather be t/t^ with fi slightly less than one (4l . 

The time evolution is believed to be caused by the growth of coherent spatial domains. 
Great importance is ascribed to the size of these domains, the coherence length %(t w ), which 
is accessible to experiments through estimates of Zeeman energies |5 1. The time evolution 
of i§(fw) plays a crucial role in the droplets theory of SG nonequilibrium isothermal dy- 
namics [6]. Perhaps unsurprisingly, it also plays a central role in yet incipient attempts to 
rationalize memory and rejuvenation experiments (see [7,89, 10 1 and references therein), 
where the experimentalist probes the glassy state by playing with the working temperature. 

Even for the simplest direct quench experiment, there is some polemics regarding the 
growth law of E, (r w ) : some theories advocate a logarithmic growth 1 6 1 , while a power law de- 
scribes numerical simulations 1 1 1 ] or experiments [ 5 ] better (a somewhat intermediate scal- 
ing has been proposed by the Saclay group 1121 and found useful in experimental work (8), 
see also Sect.[6]below). Nevertheless, two facts are firmly established: (i) the lower T is, the 
slower the growth of (r w ) and (ii) ~ 100 lattice spacings, even for T ~ 7" c and f w as large 
as 10 4 s 0. Hence, the study of SG in thermal equilibrium seems confined to nanometric 
samples, or to numerical simulations. 

There is clear evidence, both experimental [13] and theoretical 11411151 , for a thermody- 
namic origin of this sluggish dynamics. A SG phase appears below the critical temperature, 
T c . Several theories propose mutually contradicting scenarios for the equilibrium SG phase: 
the droplets [ 16 1, replica symmetry breaking (RSB) [ 17], and the intermediate Trivial-Non- 
Trivial (TNT) picture [18|. Even if this equilibrium phase is experimentally unreachable 
(at least in human time scales), we now know 1191 that it is nevertheless relevant to the 
nonequilibrium dynamics probed by experiments. 

Droplets expects two equilibrium states related by global spin reversal. The SG order 
parameter, the spin overlap q (precise definitions are given below in Sect. E), takes only 
two values q = ±qEA ■ In the RSB scenario an infinite number of pure states influence the 
dynamics 117112011211 . so that all — <7ea<<7<<7ea are reachable. TNT 1181 describes the SG 
phase similarly to an antiferromagnet with random boundary conditions: even if q behaves 
as for RSB systems, TNT agrees with droplets in the vanishing surface to volume ratio of the 
largest thermally activated spin domains (i.e. the link-overlap defined below takes a single 
value). 
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Droplets' isothermal aging [6| is that of a disguised ferromagnetQ Indeed, superuni- 
versality, the emerging picture of isothermal aging, has been found useful for the study of 
basically all coarsening systems. For T <T C the growing domains are compact geometrical 
objects. Even if the surface of these domains might be fractal, their surface to volume ratio 
vanishes as (f w ) diverges, see Eq. d 1 3 b below. Inside them, the spin overlap coherently takes 
one of its possible equilibrium values q = ±<7ea- Time dependencies are entirely encoded 
in the growth law of these domains, since correlation functions (in principle depending on 
time and distance, r) are universal functions of r/£j(f w ). 

We are not aware of any investigation of the dynamical consequences of the TNT picture. 
Nevertheless, the antiferromagnet analogy suggests that TNT systems will show coarsening 
behavior. 

As for the RSB scenario, equilibrium states with a vanishing order parameter q = do 
exist. Hence, the nonequilibrium dynamics starts, and remains forever, with a vanishing 
order parameter. Furthermore, the replicon, a Goldstone mode analogous to magnons in 
Heisenberg ferromagnets, is present for all T <T C 1 22 ]. As a consequence, the spin overlap 
is expected to vanish inside each domain in the limit of a large E, (r w ). Furthermore, q is not a 
privileged observable (overlap equivalence 1201 ): the link overlap displays equivalent Aging 
behavior. 

In order to be quantitative, these theoretical pictures of nonequilibrium dynamics need 
numerical computations for model systems. Indeed, several investigations have been carried 
out even for the simplest cooling protocol, the direct quench ll9ll 1 Oil 1 1 Il23ll24ll25 Il26ll27ll28l 
[29) ■ A major drawback of this approach, however, is the shortness of the reachable times. In- 
deed, one Monte Carlo Step (MCS) corresponds to 10~ 12 s |TJ. The experimental scale is at 
10 14 MCS (~ 100 s), while typical nonequilibrium simulations reach ~ 10~ 5 s. The problem 
has been challenging enough to compel physicists to design high-performance computers 
for SG simulations 130113 1 11321 

The situation has dramatically changed thanks to Janus [32], an FPGA computer that 
allows us to simulate the dynamics of a reasonably large SG system for eleven orders of 
magnitude, from picoseconds to tenths of a secondo Thanks to Janus, we have recently per- 
formed a study of the nonequilibrium dynamics of the Ising Spin Glass [ 331- We introduced 
novel analysis techniques that allow the computation of the coherence length in a model 
independent way. This was crucial to obtain evidence for a replicon correlator. Furthermore, 
we showed how to investigate overlap equivalence and presented evidence for it. 

In this work, we shall concentrate on the simplest protocol, the direct quench, for an 
Ising SG. We present a detailed study of dynamic heterogeneities, an aspect untouched upon 
in 1331 in spite of its relevance [27 , 28. 29 1. We show the first conclusive numerical evidence 
for a growing correlation length in the nonequilibrium dynamics, and its relationship with 
the coherence length E, (r w ) is explored. Furthermore, we compute the anomalous dimension 
for the two-time, two-site propagator (see definitions below). Due to their central role, a 
systematic way of extracting coherence (or correlation) lengths from numerical data is called 
for, but it has scarcely been investigated in the past (see however [10, 1 1 23 1). This is why 
we take here the occasion to give full details on our integral estimators 1331 (see also 1231 ). 

The layout of the rest of this paper is as follows. In Sect. [2] we define the model as 
well as the correlation functions and time sectors. We describe our simulations, which have 
been extended as compared with 1331 and discuss the difficult topic of extracting the best 



1 However, when the temperature is varied, droplets theory predicts a more complex behavior for SGs than 
for ferromagnets, due to temperature chaos. 

2 The wall-clock time needed for this computation on a cubic lattice of 80 lattice spacings is some 25 days. 
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fit parameters from extremely correlated data. Sect. [3] is devoted to the integral estimators 
of the coherence (or correlation) lengths. To the best of our knowledge, the investigation 
of this technical (but crucial) issue was started in the context of lattice field theories [34|. 
These integral estimators were instrumental to develop modern Finite Size Scaling tech- 
niques for equilibrium critical phenomena 1 35 , 36] and, therefore, to establish the existence 
of the Spin Glass phase in three dimensions [14,15,37]. In the context of nonequilibrium 
dynamics, new aspects (and opportunities) appear. In Sect. [4] we investigate the dynamic 
heterogeneities. In Sect. [5] the information gathered on length scales is used to analyze time 
correlation functions (and to extrapolate them to infinite time). We also study the thermore- 
manent magnetization. The crucial issue of the time growth of the coherence length E, (f w ) is 
considered in Sect. [6] We present our conclusions in Sect. [7] 

2 Model, correlation functions, time sectors 

2.1 Model 

The D = 3 Edwards-Anderson Hamiltonian is defined in terms of two types of degrees of 
freedom: dynamical and quenched. The dynamical ones are the Ising spins a x =H, which 
are placed on the nodes, x, of a cubic lattice of linear size L and periodic boundary condi- 
tions. The nondynamical (or quenched) ones are coupling constants assigned to the lattice 
links that join pairs of lattice nearest neighbors, J xy . In this work J xy = ± 1 with 50% prob- 
ability. Each assignment of the {J xy } will be named a sample, and, once it is fixed, it will 
never be varied [IJ. 

The interaction energy for the spins is 

Jf = - JxjO x Oy, (1) 

M 

((■■■) denotes lattice nearest neighbors). The choice J 2 = 1 fixes our energy units. We work 
in the so-called quenched approximation: any quantity (either a thermal mean value or a time 
dependent magnitude, see below) is supposed to be computed first for a given assignment of 
the couplings. Only afterwards the resulting value is averaged over the J xy , which we denote 

by— 

The spins evolve in time with Heat-Bath dynamics (see, e.g., 1381 ), which belongs to the 
universality class of physical evolution. The starting spin configuration is taken to be fully 
disordered, to mimic the experimental direct quench protocol. 

Note that the Hamiltonian (TJ has a global Z2 symmetry (if all spins are simultaneously 
reversed a x — > — O x the energy is unchanged). This symmetry gets spontaneously broken in 
three dimensions upon lowering the temperature at the SG transition at T c = 1.109(10) [39]0 

Finally, let us recall that the average over the coupling constants induces a non dynamical 
gauge symmetry 1401 . Let us choose a random sign per site e x = ±1 . Hence, the energy 

3 In principle, one should average many thermal histories for each sample and then compute the average 
over the disorder. In practice, errors can be largely reduced if one generates two histories per sample but 
generates more samples, since both sources of error are statistically independent. The claims on lack of self- 
averageness made in Sect. l2.2l depend critically on this choice. In fact, were one to simulate an infinite number 
of thermal trajectories per sample, the spin glass susceptibility would be self-averaging. 

4 The alert reader will recall that the gauge symmetry of Eq. {2} forbids a spontaneous magnetization. 
Indeed, see Sect. 12. 21 one needs to introduce an independently evolving copy of the spin configuration which 
promotes the symmetry to Z2 x Z2 . This is the symmetry which is actually spontaneously broken. 
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is invariant under the transformation 

s x * £x$xi 

i > F F i W 

Jxy t &x&yJx.y 

Since the gauge transformed couplings e x £yJ x .y are just as probable as the original ones, the 
quenched mean value of an arbitrary function of the spins 0({s x }) is identical to that of its 
gauge-average Y.{e x =±\} 0({e x s x }) /2 l , which typically is an uninteresting constant value. 
Constructing non trivial gauge-invariant observables is the subject of the next subsection. 



2.2 Observables 

A standard way of forming operators that are gauge-invariant under (O is to consider real 
replicas. These are two statistically independent systems, {o x ^} and {cj 2 '}, evolving in 
time with the very same set of couplings. Their (obviously gauge-invariant) overlap field at 
time f w is 

^(r w ) = ai 1) (f w )ai 2) (f w ). (3) 

A slight modification consists in using just one of the real replicas, say {o x ^}, but 
considering times £ w and t + r w 

c x (t , t w ) = cri 1 } (r + f w ) crj 1 ) (r w ) . (4) 

In many of the quantities defined below using c x (t,t w ), one may obviously gain statistics by 
averaging over the two real replicas. We have done so whenever it was possible, but this will 
not be explicitly indicated. We consider three types of quantities: 

1. Single-time global quantities: 

- Time-dependent energy density (N = iP is the number of spins in the lattice) 

e(«w) =-t;£ hyOx(K) e y M ■ (5) 

Recall that (x,y) indicates summation restricted to lattice nearest neighbors. 

- The spin glass susceptibility #sg(/w) is defined in terms of the SG order parameter 



Of course, the quenched mean value q(t„) vanishes in the nonequilibrium regime 
where the system size is much larger than the coherence length £(f w ) • The suscep- 
tibility 

ZSG(fw) = Nq 2 (t w ) , (7) 

steadily grows with the size of the coherent domains. Note that fluctuation-dissipation 
relations imply that ^sg ls basically the nonlinear magnetic susceptibility. 



2. Two-times global correlation functions: 
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- Spin-spin correlations: 

C(?,*w) = t;£ c x (t,ty,). (8) 

X 

The function C(r,r w ) carries many meanings: 

(a) If the first argument f w is held fixed, and C(r,r w ) is studied as t grows, it is 
just the thermoremanent magnetization. Indeed, because of the symmetry (O 
the uniform configuration that would have been enforced by holding the spin 
glass in a strong external magnetic field can be gauged to the spin configuration 
found at time r w after a random start. 

(b) On the other hand, in the pseudoequilibrium regime t -C fw> the (real part of 
the) magnetic susceptibility at frequency CO = 2k /T is given by the fluctuation- 
dissipation formula (l — C(f ,t w )) /T . 

(c) Another point we shall be concerned with is the computation of the SG order 
parameter. It may be defined from the translationally invariant time sectoi0 

C(()=limC(Mw), (9) 

as 

4EA = limCoo(f). (10) 

The computation of i^ea is notoriously difficult 1261 . Note that other authors 1271 
subtract ^ea from C<x,(t) in such a way that it tends to zero for large t. 

- The link correlation function 



Qink(Mw) = 52 c i( f . f w)cj,(f,f w ) , (11) 

carries information on interfaces. Indeed, consider a coherent spin-flip in a domain 
half of the system size. This will induce a dramatic change in C(f,f w ). On the other 
hand, the change in Cjj n k will be concentrated at the lattice links that are cut by 
the surface of the flipped domain. If the geometry of this flipped region is that of a 
compact object with a vanishing surface to volume ratio, Ci; n k will remain basically 
unchanged. 



3. Space dependent correlation functions: 
- Single time correlation function: 



C 4 (r,f w ) = T;£fe(fwk*+r(fw) • (12) 

V X 

The long distance decay of C4(r,f w ) defines the coherence length: 

C 4 (r,f w )~i/(r/!(f w )). (13) 



5 By analogy, one may define a \i time sector by the limit C£ (s)= lim fw _,„» C(s(w , t„ ) . The translationally 
invariant sector is just /J = 0. In the range < s < °° the correlation function varies in ^ea < cif -0 ' (s) < 1 . 
Full Aging 1 3 1 would imply that when < s < <*>, ci? -1 ' (.v) goes from <?ea to 0. At the present moment, it is 
unclear whether the full range of variation of the correlation function, < C < 1 , may be covered with just 
two time sectors. 
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The exponent a is relevant, because C4 at distances <§(f w ) tends to zero as <§(f w ) _fl . 
For coarsening systems, because a = 0, the order parameter does not vanish inside 
a domain. For the Ising SG in three dimensions the exponent was found to be a « 
0.4 13 3 114 1 II (see Sect. l3.2l for details). The long distance damping function / seems 
to decay faster than exponentially, f(x) = exp[— xfi] with j3 ~ 1.5 1 1 Oil 1 1 H - Note as 
well that, at the critical point, a is related to the anomalous dimension, the latest 
estimate being a(T c ) = 1 + rj = 0.625(10) 1391 . The physical origin for a nonzero 
a below T c is in the replicon mode. In fact, it was conjectured that for all T <T C , 
a(T) = a(T c ) /2 1221 . In 1331 we found values not far from this prediction. Note that 
the exponent a is discontinuous at T c |45|. 
- The two-time spatial correlation function (see I27II28I ) 



C 2 +2(r,Mw) = TT^[ c '( f ' ? w)c I+r (f,fw) - C 2 (f,f w )] , (14) 



(one could also subtract C(f,f w ) , but due to the self-averaging character of C this 
leads to the same thermodynamic limit). This correlation function is rather natural 
for the structural glasses problem, see for instance |42|, where an adequate order 
parameter is unknown. 

There is a simple probabilistic interpretation of C2+2- Let us call a defect a site 
where c x (r,r w ) = —1 and let «(f,f w ) be the density of these defects. We trivially 
have C(/,/ w ) = 1 — 2w(?,? w ). The conditional probability of having a defect at site 
x+r knowing that there already is a defect at site x is n (t , ? w ) g (r ) , where the defects ' 
pair-correlation function is g(r). Hence, C2+2(f,f,f w ) is just4« 2 (/,/ w )(g(r) — l). 
The long distance decay of C2+2{r,t,t w ) defines the correlation length £(f,? w jj 

C 2+ 2(r,f,fw)~^£(r/C(Mw))- (15) 

Basically nothing is known on exponent b nor on the long distance damping func- 
tion g. In 1 27 28 1, this decay was fitted with b = and g(x) = e~ x , but the smallness 
of the found correlation lengths £(f,y) < 2 for t + f w < 1.3 x 10 8 I27II28I , does 
not permit strong claims. In the structural glasses context [42], one tries to interpret 
£(?,f w ) as a coherence length such as (r w ), rather than as correlation length. As we 
shall empirically show, this might be very reasonable in the limit t ^> r w . 
In an RSB framework, the relaxation within a single state corresponds to the range 
^ea < C(?,? w ) < 1 (the further decay of C(?,? w ) corresponds to the exploration 
of new states). This regime is quite naturally identified with the condition that 
£(?,f w ) <C <§(?w)- In fact, <?ea yields the (correlated) percolation threshold for de- 
fects. 

Sometimes we will find it useful to change variables from t to C. This is always 
feasible, because C(t,t„) is a monotonically decreasing function of t for fixed r w . 
The accuracy of our numerical data allows this change of variable without difficulty 
(we have used a cubic spline, since the function C(t,t w ) was sampled at a selected 
set of times). 

6 The difference between coherence and correlation length is a subtle one. In this work, we shall reserve 
the name 'coherence length', which is computed from a non-connected correlation function, Eq. )12t . for the 
typical size of the coherent domains. The correlation length, which is computed from a connected correlation 
function, Eq. U4t . refers to the characteristic length for defect correlation. In particular, the coherence length 
diverges when f w — > °°, while the correlation length may or may not diverge in that limit. 
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Finally, note that C2+2 is the difference of two statistically correlated quantities 
(hence, the statistical error in the difference may be expected to be smaller than 
that for each of the two terms). This can be adequately taken into account by means 
of a jackknife procedure (see, e.g., [38 1). 

All the quantities defined so far are self-averaging (i.e., their relative errors for a fixed 
number of samples decrease as A^ 1 / 2 ), with the notorious exception of Xsg(^w)- This fact 
provides justification for the standard strategy in nonequilibrium studies (both numerical 
and experimental!) of averaging results over very few samples. 

Self-averageness stems from the fact that the computed/measured quantities are averages 
of local observables taken over the full system (which provides a number of statistically 
independent summands of the order of L D /<§ D (f w ). The exception, Xsg(/w) Q, is actually 
non local as it is the integral over the whole system of C4(r,/ W ), Indeed, the central limit 
theorem suggests that the probability distribution function of q{t w ) should tend to a Gaussian 
when L — > °°. Hence, the variance for Xsg(?w) is ~ 2#| G (f w ) in the limit of a large system. 

2.3 Simulation details 

The Janus computer [ 32] can be programmed for the simulation of the single spin flip Heat 
Bath dynamics up to a very large number of lattice sweeps (units of Monte Carlo time) for 
systems of linear sizes up to L ~ 1000 

We have spent the most effort in simulating the dynamics of the model described by 
Eq. (TJ in the direct quench protocol described in the Introduction, for several runs of about 
a hundred samples of linear size L = 80 and up to 10 11 Monte Carlo steps (we recall that 
a single step corresponds to roughly 1 picosecond of time in the real world). Details of our 
simulations are given in Table [TJ We extend here the analysis of the simulations reported 
on 1331 . but additional simulations have also been carried out. Most notably, we simulated 
768 new samples of size L = 80 at T = 0.7 up to 10 10 in Monte Carlo time, which have been 
useful to improve and test the statistical accuracy in some aspect of our analysis. 

We wrote to disk the spin configurations at all times of the form [2 ,//4 ] + [2 J / 4 ], with 
integer and j (the square brackets stand for the integer part). Hence, our t and r w are of the 
form [2'' 4 ]. Nevertheless, we computed C2+2 only for powers of two, due to the increased 
computational effort. 

A final note on the time span of our runs. Much to our surprise, we found in [33 ] that, 
even for our very large systems, finite size effects in the coherence length can be resolved 
with our statistical accuracy. In this work, we have restricted ourselves to the time window 
that is not affected by them. The single exception will be in the analysis of energy relaxation, 
Sect. 15.21 where this range is too short. Nevertheless, we have explicitly checked that the 
energy suffers from smaller finite size effects than the coherence length. 

2.4 Fits for extremely correlated data 

Computing the best fit parameters and estimating errors from extremely correlated data sets 
presents a common, and still not satisfactorily solved, difficulty in many numerical studies. 
For instance, our study of C(r,r w ) requires considering approximately 10 4 random variables 



7 The overall parallel update rate for L = 80 systems of the whole Janus — 256 nodes — is 78 femtoseconds 
per spin. 
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L 


T 


MC steps 


N s 


80 


0.6 


10" 


96 


80 


0.7 


10" 


63 


80 


0.8 


10" 


96 


80 


0.9 


2.8 x 10 10 


32 


80 


1.1 


4.2 x 10 9 


32 


80 


1.15 


2.8 x 10 10 


32 


80 


0.7 


10 10 


768 


40 


0.8 


2.2 x 10 s 


2218 


Table 1 


Parameters of our simulations. 


The overall wall-clock time needed 


was less than six weeks. We 



highlight with boldface the simulations performed after completion of [33|. Recall that we take the critical 
temperature from [39], T c = 1.109(10). The full analysis of spin configurations was performed offline. 



extracted from a set of only 63-768 samples. The standard approach, computing the covari- 
ance matrix and inverting it, fails because this matrix is necessarily singular^ In this paper 
we shall follow an empirical procedure. We shall consider only the diagonal part of the co- 
variance matrix in order to minimize x 2 when performing fits. Unless otherwise indicated, 
we shall always use this diagonal x 2 in the rest of the paper. Yet, in order to take correlations 
into account we shall perform this procedure for each jackknife block and later on compute 
error estimates from their fluctuations. As we have run simulations with both 63 (from | 33]) 
and 768 samples for T = 0.7, we are in a position to test this method by comparing the 
results obtained and the ones to be expected for 63 samples (see Sect. l3.2t . 

The main drawback of this approach is that the standard x 2 test of fit likelihood cannot 
be applied blindly. Of course, were the exact fitting function known, the average value of 
diagonal x 2 should be 1 per degree of freedom. Yet, since the obtained fitting function may 
coherently fluctuate with the numerical data, we shall encounter anomalously low values of 
diagonal x 2 - We examine this problem in Sect. 13.21 empirically finding that x 2 behaves as 
if there were many fewer degrees of freedom than what one would expect. 

3 Integral estimators of characteristic length scales 

The need to estimate characteristic length scales, such as ^(f w ) or £(r,r w ) is a recurrent 
theme in lattice gauge theory and statistical mechanics. The more straightforward method is 
to consider a particular functional form for the long distance damping function in Eqs. d 1 3 b 
or d 1 5 b - One of the problems with this approach, already identified in the study of equi- 
librium critical phenomena, is that it is extremely difficult to extract from numerical data 
simultaneously the length scale and the exponent for the algebraic decay. Note that quite 
often computing the exponent is as important as extracting the length scale to draw physical 
conclusions. The situation worsens if the functional form is only an educated guess, which is 
precisely our case. Furthermore, numerical data for the correlation function at different lat- 
tice sites suffer from dramatic statistical correlations, which complicates fitting procedures. 

A different approach, the use of integral estimators, has been known since the 1980s [ 34 1; 
but only in the mid 1990s (see e.g. 1 35 .36]) it was realized that it provided an enormous 

8 Consider computing the covariance matrix for a set of No random variables from JV S < No samples. If 
the N s x No numbers are disposed on a rectangular matrix, it is clear that the last No — N s rows (reordering the 
rows if necessary) are a linear combination of the others. In other words, the results are indistinguishable from 
the ill-conditioned situation where these last No — N s random variables are the very same linear combination 
of the first ones. Once this is realized it is trivial to show that the range of the size No x No covariance matrix 
is at most N s . 
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simplification. The use of integral estimators for the length scale enables determinations of 
exponents such as a in Eq. d 1 3 b - which are completely independent from the functional form 
of the long distance damping. The only place left for systematic errors is in finite size effects 
or in scaling corrections (when the considered range for the variation of length scales such 
as £, (f w ) is too small). As for the determination of the length scale itself, integral estimators 
are guaranteed to produce numbers that scale as the inaccessible true | (f w ), provided that it 
is large enough. 

The fact that the correlation functions that will be considered here, C\ and C2+2, are 
self-averaging in a nonequilibrium context provides an impressive error reduction, which is 
not accessible for equilibrium studies. 

Our chosen example to explain the method will be that of C4 and the determination of 
the coherence length and of the exponent a. 



3.1 The coherence length 



Cooper et al. | 34 1 suggested the second moment determination of the characteristic length, 

n 1/2 



^ (2) ('w): 



1 



V^sin tt/L 



C 4 (0,? w ) 

C*4 (^min? A* 



1 



(16) 



where C4(k,t w ) is the Fourier transform of C4.(r,f w ), and ft m ; n is the minimal non- vanishing 
wave vector allowed by boundary conditions (ft m j n = (2jt/L,0,0) or permutations). Notice 
that Xsc{tw) = £4(0, t w ). As can be readily seen, in the thermodynamic limit this is equiva- 
lent to 



JdOrC 4 (r,f w ; 



(17) 



The denominator in this equation is just the SG susceptibility I0 which, as we said in 
Sect. 12.21 does not self-average (and neither does the numerator). Because of this, if one 
were to follow this method, a very large number of samples would be needed. 

We would like to have a better statistically behaved definition of £; (f w ). In order to get 
it, we start by considering the integral^ 



40w 



L/2 



drr*C 4 (r,f w 



(18) 



As we are going to work in the thermodynamical limit, we are interested in the regime 
L>^ (f w ), so we can safely reduce the upper integration limit from °° to L/2. 

With this notation and assuming rotational invariance, the second moment coherence 
length is just 



Ip+\ (tyj) 
fo-l(fw) 



(19) 



9 In what follows C 4 (r,t w ) stands for C 4 (r,^), with r = (r,0,0) and permutations. As we shall see in 
Sect. |3.3| using an average over spherical shells does not achieve a significant reduction of statistical errors 
in our chosen estimator for the coherence length. 
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Fig. 1 The spatial autocorrelation of the overlap field for f w = 2 and three subcritical temperatures, as 
computed in our L = 80 lattice. On the left panel we show r 2 Ct(r,f w ), recall while on the right one we 
show r 4 C4(r,f w ) (mind the different scales). While the signal to noise ratio of both quantities falls equally 
rapidly, the problem is less severe for the computation of I2 (t w ) since the maximum there is not in the noise 
dominated region. The curve for T = 0.7 is the average of 768 samples, while those of T = 0.6 and T = 0.8 
are computed from 96 samples. 



We also recall that in (23l it was proposed to identify | (r w ) with /o(/w)0 but this would 
only be appropriate for a = 0. For a correlation function following the scaling law d!3l l. one 
can use a more general definition, because 4(f w ) x [£, (t w )] k+i ~ a : 

| M+1 ( fw ) = ^^oc|(f w ). (20) 

Definitions such as d!6t and d20l i suffer from systematic errors because equation d 1 3 b is only 
an asymptotic formula for large values of r. Therefore, the systematic errors in these defi- 
nitions can be reduced by considering a large value of k (since the r k factor would suppress 
the deviations at short distances). However, there is also the issue of statistical errors to con- 
sider. As we can see in Fig.[T] a large power of r pushes the maximum of r^C^r, t w ) into the 
region where r 2> £, (f w ) and the signal to noise ratio of the correlation function is extremely 
low. Because of this, a compromise in the choice of k is needed. Our preferred option is 
ll,2(«w)- 

Even though our use of §1,2 (?w) instead of ' 2 ' (f w ) already mitigates the statistical prob- 
lems in the integration of G4.(r,J w ), we can still improve the computation. As can be plainly 
seen in Fig. [TJ r 2 C4(r,t w ) starts having very large fluctuations only at its tails, where the 
contribution to the integral is minimal. To take advantage of this fact, we are going to use 
a self-consistent integration cutoff (a method applied before in the study of correlated time 
series |43]). We only integrate our datJH for C4(r,r w ) up to the point where this function 
first becomes less than three times its own statistical error. Of course, while this method 
provides a great reduction in statistical errors, it does introduce a systematic one. To avoid 
it, we estimate the small contribution of the tail with a fit to 

Q(r,«w) = ^4 exp [-(r/| fit (f w )) L5 ] ■ (21) 

10 The authors of 1231 carefully discussed the interplay of the integration limits in )18t and the boundary 
conditions. Since we restrict ourselves to L 3> £(f w ), this is immaterial to us. 

As our (somewhat arbitrary, yet irrelevant) choice of quadrature method we have chosen to interpolate 
the data with a cubic spline, whose integral can be exactly computed. 
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Fig. 2 Left: Result of computing £i ? in two different ways for our 96 samples at T = 0.6. In the orange 
curve (color online) we stop the integration at the cutoff point where relative error of C4 is greater than one 
third. In the blue curve (color online) we estimate the contribution of the tail from that point on extrapolating 
with a fit to <2 1 1 . The difference is small, but with the second method the power law behavior of <^i.2('w) lasts 
longer. Right: Same plot for our 768 samples at T = 0.7. With the increased statistics this extrapolation is 
not as important and both curves are compatible for the whole simulation. With the 63 samples of [33], the 
tail contribution is as significant as in the left panel. 



Notice that this is just the scaling function i 13k using f(x) = exp[— x 15 ] as our damping 
function and a = 0.4. Of course, while this fit is used to estimate the contribution of the 
interval [r cuto ff,L/2], we actually perform the fitting for 3 < r < min{15,r cuto ff}, where the 
signal is still good. This last step is important for large r w (see Fig. [2}. 




Fig. 3 The SG susceptibility Xsg(<w) for our 768 samples at T = 0.7 computed from q 2 and from the integral 
h - ZsG(fw) — Nq 2 (t w ) = 47Ch(t w ). The main difference between the two determinations is that the second 
one has been computed with a self-consistent cutoff. As we can see, even though both curves are compatible, 
the integral one is much more precise. The inset details the upper right corner. 



As a consistency check of this method and as a demonstration of its enhanced precision 
we can consider the SG susceptibility I0. This observable, Xsg(/w) = Nq 2 (t w ), coincides 
with 4^/2 (? w ) in the presence of rotational invariance. We have plotted both expressions as a 
function of time in Fig. [3] The only systematic discrepancy between the two is at short times, 
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Fig. 4 Comparison of our integral estimators ijcu and §12, Eq. )20t . with the second moment estimate i;' 2 ' 
and the result § flt of a fit to i2l"l , with a = 0.4. All the curves become parallel at large r w , but the integral 
estimators have much smaller errors. All curves are for our 768 samples at T = 0.7. 



With tails ■ 63 samples 

Without tails 768 samples 



Fig. 5 Left: Ratio of the errors in £1 2(Av) for the simulations at T = 0.7 with the 63 samples of [33 1 and 
those of our simulations with 768 new samples (see text for discussion). The extrapolation to include the tails 
in the integrals is immaterial for this ratio. The horizontal line is ^768/63 ~ 3.46. Right: Cutoff of the 7j. 
integrals as a function of time for both simulations. 



when the system cannot be considered rotationally invariant (see Sect. 13 - 3b - However, the 
integral determination 4nh (fw) is much more precise for the whole span of our simulation. 

As a second check, we have plotted in Fig.|4]the integral estimators |o,i ( f w) and (^^(fw), 
together with the traditional second moment estimate ^' 2 ' and the result of a fit to d2 1 b - As 
we can see, all determinations are indeed proportional, but the integral estimators are much 
more precise. 

Finally, we address the issue of our error estimates by comparing our data at T = 0.7 
for the 63 samples of [ 33 1 with our new simulations with 768 samples. We have explicitly 
checked that the errors in Cn(r,t m ), computed with the jackknife method, scale as the inverse 
of the square root of the number of samples, within errors (for Gaussian distributed data, the 
relative statistical error in the error estimate is ~ l/^/2iVsampies)- The fact that C4 verifies 
this basic expectation is a demonstration that large deviations, which would be missed for 
a small number of samples, do not appear, even for r w as large as 10 10 . On the other hand, 
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Fig. 6 Difference between the coherence length ^i.2(?w) computed with the 63 samples of ['33] and with the 
768 new samples (the errors are the quadratic sum of those for each simulation). Both curves are compatible 
in the whole time range. Mind the dramatic statistical correlation in the sign of this difference. 

the behavior of statistical errors in integrals with a dynamically fixed cutoff is not that sim- 
ple 1431 . In our case, the ratio of the errors in 

|i,2(fw) (Fig- El left) is around 30% below 
the Gaussian expectation. We have checked that a similar effect arises in the computation of 
h(tw) and that the effect of the tails is immaterial. In fact, this deviation is entirely due to 
the difference in the dynamical cutoffs of both simulations, Fig.[3J right. Whenever the cut- 
off coincides, the error ratio is in the expected region around a/768/63 « 3.46. The overall 
consistency of our error determinations is demonstrated by Fig. [6] 



3.2 The algebraic prefactor 

One of our main goals is to provide a precise estimate of the exponent a for the algebraic 
part of C4, see equation d 1 3 b - Equilibrium methods [36] are not well suited to a nonequilib- 
rium study in the thermodynamic limit. Instead, we introduce here the method used, but not 
explained, in 1331 . 

The starting point is the realization that 1\ °c <§ 2 2°, which would indicate that a can be 
obtained from a power law fit of I\ as a function of the coherence length. Furthermore, the 
large statistical correlation between I\ and |i 2 can be used to reduce the statistical errors 
in a. However, such a fit would be quite problematic, as we would have errors on both 
coordinates (the correlation of the data already poses a nontrivial problem with errors in just 
one coordinate, see Sect. 12. 4\ . Instead, we fit separately I\ (f w ) and §1 2(*w) to power laws in 
the waiting time. This way, if 

/!(r w )=A4, $(t w )=Bt l J z , (22) 

we have a = 2 — cz. This relation holds for each jackknife block, which lets us take full 
advantage of the correlations. 

Following this method, we obtain the results in Table [2] In [33 ] (first four rows of Ta- 
ble[2]l we quoted the results for a fitting range of ^ 2 € [3, 10], which is perfectly adequate 
for T = 0.6,0.8 and 1.1. As we can see, x 2 appears to be a bit too large for T = 0.7, but if 
we narrow the fitting range it becomes reasonable and a does not change. 
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T 


^ samples 


[^?min i ™>inaxj 


z 


a 




2/j/d.o.f. 


0.6 


96 


[3,10] 


14.06(25) 


0.359(13) 


41.7/82 


49.0/82 


0.7 


63 


[3,10] 
[3.5,10] 


11.84(22) 
12.03(27) 


0.355(15) 
0.355(17) 


82.7/81 
52.7/71 


131/81 
75.5/81 


0.8 


96 


[3,10] 


9.42(15) 


0.442(11) 


17.1/63 


12.2/63 


1.1 


32 


[3,10] 


6.86(16) 


0.585(12) 


18.7/46 


26.1/46 


0.7 


768 


[3,10] 
[3.5,10] 

[4. 10] 
[4.5, 10] 


11.45(10) 
11.56(13) 
11.64(15) 
11.75(20) 


0.395(8) 
0.397(10) 
0.397(12) 
0.394(14) 


86.9/76 
46.6/66 
40.1/58 
29.6/50 


269/76 
101/66 
60.4/58 
35.8/50 



Table 2 Value of the dynamic exponent z and the algebraic prefactor a for several temperatures. The fitting 
range §12 6 [3,10], which worked for the smaller number of samples we had in 1331 . does not give good 
fits for /[ with our enlarged statistics at T — 0.7 (the fits for the coherence length itself are still good). 
Nevertheless, if we increase ^ m i n to get reasonable values of /d.o.f. we see that the estimate of a does not 
change. 
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Fig. 7 Left: Probability density function p(a) of the estimate of exponent a, Eq. 4 13K as obtained from a 
set of 63 samples. The dots with horizontal error bars are, from top to bottom: our best estimate with 768 
samples, the value with the 63 samples of [33], and the mean and standard deviation of p(fl). The continuous 
line is a Gaussian distribution with the same mean and variance as p(a). Center: As in left panel, for the 
jackknife errors Aa. The vertical line marks the standard deviation of p(a). Right: Histogram of the ^ 2 /d.o.f. 
parameter for the xi fit. In all three panels the fitting range to obtain a was taken as t, £ [3, 10]. 



As this is a very important magnitude, in this work we have increased the number of 
samples at T = 0.7 by a factor of 12 with respect to 1331 . (from 63 to 768 samples, with 
extra simulations that stop at f w = 10 10 , rather than 10 11 ). This not only allows us to provide 
a better estimate of a but we are now also able to check the soundness of the statistical 
procedure. 

The first difficulty is that, with the corresponding reduction in statistical errors, the orig- 
inal fitting window no longer provides reasonable values of our diagonal % 2 estimator. In- 
stead, we have pushed the lower limit to | > 4 (see Table|2]for details). The new value of 
a(T = 0.7) is 

a(T = 0.7) = 0.397(12). (23) 



In accordance to the analysis of Fig. [5] the error has not decreased the factor 1^/768/63 
with the increase in statistics. One could say that the dynamic cutoff procedure has traded 
statistical uncertainty for a reduction in systematic errors. Another contributing factor to the 
large statistical error is the raising of the minimal coherence length included in the fit. 
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To understand whether the discrepancy between the new estimate of a and the one in 1331 
is due to a systematic effect or to a large fluctuation, we can use a Monte Carlo method. The 
probability distribution function of the estimates of a, as computed with 63 samples, can 
be obtained easily from our set of 768 samples. One randomly picks sets of 63 different 
samples and determines a and its jackknife error Aa for each of these sets (mind that there 
are ( 7 g 6 3 8 ) w 2.2 x 10 93 possible combinations). We have done this 10 000 times and computed 
normalized histograms of both quantities (Fig.UJl. Clearly enough, the estimate in 1 33 ] was a 
fluctuation of size 2.2 standard deviations, large but not unbelievably so. On the other hand 
we see that the jackknife method tends to slightly underestimate Aa for 63 samples (Fig. [7] 
center). Note as well that there seems to be a small bias (smaller than the error) on the 
estimate of a with only 63 samples (Fig. [JJ left, compare the histogram with the uppermost 
horizontal point). There are two possible reasons for this. One is that a is obtained from raw 
data through a nonlinear operation. The other is that the larger the number of samples, the 
smaller the cutoff effects in the computation of 4(f w ). 

It is amusing to compute as well the probability density of ^ 2 /d.o.f. for fits with 63 
samples. We show in Fig.[7J right, that %^ for the fit of the coherence length can be much 
larger than what one would naively expect for a fit with ~ 80 degrees of freedom. 



3.3 Isotropy of C4(r,t w ) 

In the previous section we neglected the issue of the isotropy of C4(r,f w ). At all times we 
worked with radial functions C4(r,f w ), obtained by averaging the correlation at distance r 
along the three axes. In doing this, we ignored most of the N points that C4(r,? w ) has for 
a given t w . The main motivation for doing this is avoiding the computation of the whole 
correlation function, a task which in a naive implementation is ff(N 2 ). Of course, due to the 
Wiener- Khinchin theorem, we can reduce this computation to the evaluation of two Fourier 

y 




Fig. 8 Level curves C4(r,f w ) = c for c = 0.3 (dashed lines) and c = 0.1 (solid lines) at T = 0.6 (dotted lines 
are circles, for visual reference). We have restricted ourselves to the z = plane, for clarity. The innermost 
curve corresponds in both cases to t w = 4 and the succeeding ones correspond to geometrically growing 
times (/ w = 4 x 16')- As we can see, the deviations from isotropy are mainly due to lattice discretization (i.e., 
functions of r), even though there is also a small dependence on time for curves of similar radii. The errors 
are smaller than the thickness of the lines (the interpolation to draw these continuous curves from our lattice 
data was performed with Mathematica TM). 
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Fig. 9 Left: Comparison of Q 2 {r,t w = 2 20 ) with r 2 C 4 (r,( w = 2 20 ) for T = 0.6. The first quantity is obtained 
by averaging over spherical shells, while the second one considers only correlations along the axes. The 
behavior of Qi is better at the tails, but does not imply any gain in practice, as both functions are equally well 
behaved up to the cutoff point. Right: The coherence length computed with the whole correlation functions, 



-JL» 



, and with correlations along the axes, £12. Both estimates coincide for large times. 



transforms which, using an implementation of the FFT algorithm 1441 . is an ff(NlogN) 
task. 

We shall examine in this section whether the complete correlation functions are isotropic 
and whether we can take advantage of them to reduce the errors in our determination of the 
coherence length. The first question is answered by Fig. [8J where we compute the level 
curves C4.(r,f w ) = c for several values of c and t w . As we can see, isotropy is recovered at 
quite small distances (remember we are only concerned with | > 3). 

In order to use our integral method for a three-dimensional Ci(r,r w ) we must first aver- 
age it over spherical shells. We do this by defining the functions Qk(n,t w ), 

£ \r\ k C 4 {r,t w ) 
Qk(n,t w ) = ^ j . (24) 

|r|e[«.«+l) 

Notice that 2o(0,r w ) = C4(0,f w ) = 1 and that the division by the number of points is needed 
to average over the spherical shell. Now we can use Qk{r,t w ) in the same way we used 
r*C4(r,f w ) in the previous section. The resulting coherence length |i 2 (Av) would be ex- 
pected to coincide with §1,2 (*w) m me large f w limit, but have much smaller errors due to the 
large increase in statistics. As we can see from Fig. |9] however, the correlation among the 
points is so great that the gain in precision is insignificant. 

We can conclude from this section that the usual approximation of considering C4(r,t w ) 
isotropic is a well founded one and that it is safe, and statistically almost costless, to restrict 
ourselves to correlations along the axes, as we did in B3I and in the previous section. Nev- 
ertheless, the computation of the whole C^/Vw) could be rewarding if one were to study 4, 
with k>2. 



4 Characteristic length scales for dynamical heterogeneities 

The concept of heterogeneous dynamics has been recently borrowed from structural glasses 
to describe non-local aging in Spin Glasses | 27 , 28 1. A coarse-grained spin correlation func- 
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Fig. 10 Left: The two-time spatial correlation function C2+2(r,f,?w) (see Eq. 1141 ) as a function of r, for a 
high value of t w = 2 35 , and several values of t. Right: C2+2 (r, t, f w ) as a function of t and several values of f w , 
fixing r = 10. Both figures are for our 63 sample simulation at T = 0.7. 
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Fig. 11 Left: Coherence length £ (C 2 ,f w ) as a function of C 2 for several values of f w at T = 0.7 (63 samples). 
Right: As figure left, for 7" = 0.8 (96 samples). 



tion may be defined for which the microscopic average is not performed on the overall vol- 
ume but on cells of size £ D . Non-uniform aging of the system results in spatial fluctuations 
of the coarse-grained correlation functions, that can be used to define a two-time dependent 
correlation length. In fact, at fixed system size and temperature one may study the two-time 
dependent distribution of the coarse-grained correlation function values or, equivalently, the 
distribution as a function of i w and of the global spin-spin correlation function C(f,f w ), 
Eq. ([SJ: if the coarse-graining size £ D is larger than the correlation length, then one should 
observe Gaussian statistics, while strong deviations from a Gaussian distribution are present 
in case of small £ D 1281 . Such dependence of the statistics on the coarse-graining size defines 
a crossover length £(C ,f w ) interpreted as an aging (two-time) correlation length. 

It has been observed 1 28 1 that such an aging correlation length may be obtained from the 
spatial decay of the two-time two-site correlation function C2+2, Eq. ( 114b : still, the authors of 
reference | 28 1 could not measure £ values greater than two lattice spacings. In this section 
we present data from our simulations on lanus showing correlation lengths for dynamical 
heterogeneities up to order ten lattice units. We show C2+2(r,f,f w ) for f w = 2 35 and some 
values of t at temperature T = 0.7 in the left picture of Fig. [TO] (as for C4, we denote by 
C2+2(r,f,? w ) the correlations along the axes). As one can see in this figure, at large times, 
correlations grow up to several lattice spacings. 
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In what follows it is convenient to eliminate the time t dependence in favor of the global 
spin-spin correlation function C(?,? w ). For given t w and C values, we obtain easily f(f w ,C) 
because the monotonic time dependence of the (discrete measures of the) correlation func- 
tion can be smoothly interpolated by means of cubic splines. Then, for each value of r, we 
must interpolate C2+2(r,t (C ,t w ) ,t w ) . As one can see in the right picture of Fig. [TOl there is a 
sharp change in the t -derivative of C2+2, so we had to resort to a linear interpolation in order 
to avoid the strong oscillations that a spline had suffered from. Once C2+2(r,f (C,f w ),f w ) has 
been interpolated at all r values, the methods of Sect.[3]allow us to estimate the correlation 
length £(CVw). 

In the large correlation sector, i.e. g| A < C 2 , and for large t w , the correlation length 
£(C 2 ,f w ) approaches a f w independent value. On the other hand, one expects that for small 
values of C 2 (that is, C 2 < (?ea) and large f w , £(C 2 ,f w ) diverges as the coherence length 
<§(f w ) defined in Eq. d 1 3 b |29|. Such behavior is represented in Fig. [IT] in which we plot £ 
as a function of C 2 , at temperatures T = 0.7 and T = 0.8, for some values of f w . It is also 
interesting to consider the ratio R(C 2 ,t w ) = £(C 2 ,f w )/^(f w ) and study how its behavior as 
a function of C 2 changes with f w . As pointed out above, for small values of C 2 and large 
f w , we expect 7?(C 2 ,f w ) ~ const > 0. Moreover, since the coherence length <§(f w ) diverges 
for large f w , R(C 2 ,t w ) should vanish at large waiting times when C 2 > g| A . In Fig. 1121 we 
show R(C 2 ,t w ) at temperatures T = 0.6, 0.7 and 0.8. An interesting feature is the crossover 
between the two sectors C 2 < g| A and C 2 > q 2 EA . At T = 0.8 the g EA is too small (see 
Table |4] in Sect. 15.11 below) to let us observe the small C 2 behavior described above. On 
the other side, data for R(C 2 ,t w ) at T = 0.6 quickly approaches a constant value for small 
correlations. It seems that the larger f w , the fastest the convergence to a constant, determining 
in this way a crossing point. However, these data suffer from large fluctuations that do not 
allow us to make any strong speculation on the crossings among curves at large values of 
t w . Indeed, we know that R(C 2 ,t w ) should vanish for large C 2 roughly as 1 /£, (r w ), but up to 
our knowledge, there is no reason for R, as a function of C 2 , to converge faster to a constant 
when f w increases. In the bottom pictures of Fig. [12] we report the same data at T = 0.7, 
averaged on both 63 (left) and 768 samples (right). Even if simulations on larger sample 
statistics are not as long as those of the first 63 samples, the smoothing in the curves does 
not improve the crossing definition. In addition, we have not been able to find any clear 
scaling behavior of the crossing points with the waiting time. 

Since at large times and small correlations £(C 2 ,? W ) and |(f w ) on ^y differ in a constant 
factor (which is also close to 1), one may expect that the behavior of the two-time, two- 
site correlators (Eq. dl4l >) should be analogous to that of the four point correlation function 
(Eq. Ill2b). We can then probe the long distance scaling of C2+2 with the same analysis 
performed in reference 1331 (and in Sect. 13.21 ) for Czt(r,f w ), Eq. d 1 3 1 > . In particular, we can 
extract the exponent b of the algebraic prefactor in Eq. d!5l l and the dynamic exponent zr, 
assuming a power-law growth for the correlation length: 

C(C 2 ,f w )=A^. (25) 

We fix C 2 to a ^-dependent value, which is a value small enough to be below the g| A at 
all considered temperatures T = 0.6, 0.7 and 0.8 (C 2 values should be compared with q 2 EA 
estimates given in Sect. l5.lb . and large enough to have the necessary number of r w points in 
order to obtain fair fits. We also had to impose a £1.2 > 3 constraint to avoid the effects of 
lattice discretization. We did not impose any upper limit to £, whose values hardly reach 8 
lattice spacings at T = 0.8 (even less at colder temperatures), and we expect that in the time 
sector considered the constraint imposed on in reference ( 33 1 and in Sect. l3.2l would work 
as well as in that case. 
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Fig. 12 Behavior of the ratio R(C 2 ,t w ) = £(C 2 ,f w )/^ ('w) as a function of C 2 for several values of f w . Top 
left: at our lowest temperature T = 0.6. Top right: at temperature T = 0.8. Bottom left: at temperature 
T = 0.7 averaged over 63 samples. Bottom right: same plot as bottom left, averaging over 768 samples. 



T 


C 1 




b 


^/d.o.f. 


XtJA.o.1 


0.6 


0.200 
0.325 


13.4(6) 
12.8(4) 


0.43(4) 
0.55(3) 


0.01/2 
7.16/7 


0.13/2 
4.45/7 


0.7 


0.200 
0.325 


11.14(20) 
11.35(12) 


0.508(17) 
0.642(9) 


0.69/3 
8.08/7 


0.37/3 
7.70/7 


0.8 


0.100 
0.175 


9.56(17) 
10.12(13) 


0.497(13) 
0.540(10) 


3.73/5 
7.15/8 


3.37/5 
7.91/8 



Table 3 Value of the dynamic exponent and exponent b for the algebraic prefactor, for three subcritical 
temperatures. We limited the fitting window constraining the fits to £12 > 3. Data at T = 0.7 averaged over 
768 samples. 



We summarize the obtained values of b in Table [3] reporting results for the smallest 
C 2 attainable (that is, permitting reasonable fits) at all temperatures, as well as for larger 
values of C 2 , for which the number of t w points allows for quite good fits. We see that b and 
Zr are slightly different from the values of a and z for the four point correlators presented 
in Sect. 13.21 Unfortunately our data do not permit a more precise determination for these 
exponents in the deep C 2 < <?| A sector. In this respect, the determination of qEA is a crucial 
issue (see Sect. l5.lb . 
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Fig. 13 The spin-spin correlation C(r,f w ) as a function of x 2 = (£(f,f w )/<!;(f w )) for T = 0.7. Inset: Close 
up of the small x 2 region and comparison with our quadratic fits. 

5 The time correlation functions 

5.1 The stationary part of C(t,t w ) 

The naive computation of the stationary part of C(t,t w ), Coo(f) Eq. (9), suffers from an es- 
sential problem: how to know when f w is large enough. The consideration of characteristic 
length scales may simplify this problem, as the limit f w — > °° is equivalent to | (? w ) 1 — ► 00 
However, the approach to this limit will be acutely f-dependent. Hence, it is better to con- 
sider a dimensionless variable, 

= (26) 

As we saw in Sect. [4] the correlation length £(r,r w ) quickly reaches a ? w -independent limit, 
so x(t, f w ) is essentially _1 (f w ) in its natural units for each t. In fact, see Fig.Q~3]for our data 
at T = 0.7, the plot of C(f , f w ) against x 2 (r, / w ) is pretty smooth for x 2 — > 0. Furthermore, the 
curves for different r become parallel as f grows, which suggests the existence of a smooth 
scaling function, C(f,f w ) = Co»(f) +/(x 2 ). We have fitted the curves C(t,x 2 ) for each t in the 
range x 2 < 0.5 to a quadratic polynomial (see inset to Fig. 113b. 

C(t,x 2 ) =C„(t)+a l (t)x 2 + a 2 {t)x 4 . (27) 

Unlike our treatment of the exponent a (Sect. [3~5l l. here we cannot skirt the issue of errors in 
both coordinates. As the effect of both errors is similar, we have performed a least squares 
fit and an a posteriori % 2 test, Fig. [14] right. This % 2 test was diagonal in the sense that 
we did not consider correlations among different pairs of times, but we did consider the 
covariance matrix for (x 2 ,C) at the same (f,f w ). Our statement that the scaling curves in 
Fig. Q~3] become parallel is made quantitative in Fig. [T4] left, which plots the coefficient 
a\(t) in equation d27| i. 

The method described above has allowed us to compute C„{t ) with remarkable accuracy 
for t < 10 , see Fig. [15] We can now try to perform a second extrapolation, Eq. d 10b . and 



In what follows we shall always use the ^12 estimator. 
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Fig. 14 Left: Slope at x 2 = of the fitting curves {27) at T = 0.7 as a function of time. Right: Check of 
the fit quality using a % 2 estimator that takes into account correlations only at equal (f,f w ), but disregards 
correlations among different pairs of times. 
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Fig. 15 Stationary component of C(t,t w ), C*,(t), obtained with the extrapolations in Figs. l 1 3l and ll4l for all 
our subcritical temperatures. 



obtain the value of q^A for each of our simulated temperatures. In order to do this we first 
tried a power law extrapolation 



Co(f) = q E A +At B 



(28) 



This functional form yielded very good fits for T = 0.6,0.8 but the values of the exponent 
B where very small, of about B ~ —0.05. Unfortunately, the smallness of B makes the ex- 
trapolation extremely risky. Just to be on the safe side and check explicitly that q^A > 0, we 
have tried a logarithmic fitting function, 



Co(r) = <?ea + 



5 + logr 



(29) 



This Ansatz lacks any theoretical basis, but since it is slower than any power law, we expect 
it to provide a lower bound on q^A- On the numerical side, the logarithmic fit was as good 
as the power law (as determined by a x 2 test). Nevertheless, as expected, it produced values 
of qEA which were incompatible with those of equation d28b (see Table |4j. Furthermore, 
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Fitting 




>arithm 




Power law 




T 


range 


4EA 


* 2 /d.o.f 




—B x 10 2 


r7 d -°-f 




[10 2 ,10 s ] 


0.607(16) 


34.1/17 


0.730(8) 


5.7(4) 


31.2/17 


0.6 


[10 3 ,10 8 ] 


0.62(3) 


7.23/14 


0.733(14) 


5.8(7) 


7.59/14 




[10 4 ,10 8 ] 


0.62(5) 


6.25/10 


0.726(24) 


5.4(12) 


6.32/10 




[ioMo 8 ] 


0.497(10) 


23.7/17 


0.656(5) 


6.16(18) 


32.6/17 


0.7 


[10 3 , 10 s ] 


0.474(21) 


18.9/14 


0.637(11) 


5.5(3) 


18.5/14 




[10 4 ,10 8 ] 


0.49(5) 


15.0/10 


0.63(3) 


5.4(9) 


15.3/10 




[10M0 8 ] 


0.371(13) 


6.50/17 


0.568(7) 


6.56(20) 


9.39/17 


0.8 


[10 3 ,10 8 ] 


0.368(24) 


5.53/14 


0.556(12) 


6.2(4) 


4.27/14 




[10 4 ,10 8 ] 


0.40(6) 


4.31/10 


0.56(3) 


6.4(11) 


3.82/10 



Table 4 Estimate of for three subcritical temperatures, using two different extrapolating functions. For 
T = 0.6,0.8 both are very good, but at T = 0.7 (where we have better statistics) they are somewhat forced. 
This suggests that the real ^ea probably lies in between our two estimates. For the power law extrapolation, 
Eq. {28), we also quote the exponent B. Notice that this exponent is not proportional to T. 



when we tried both extrapolating methods for T = 0.7, where we had simulated many more 
samples, we found that they were somewhat forced. This leads us to conclude that the real 
asymptotic behavior of C^(t) is probably something in between equations d28l l and d29l l. We 
shall use the difference between both methods with a fitting window of / 6 [10 3 , 10 8 ] as our 
uncertainty interval, 

0.62 < q nh {T = 0.6) < 0.733, 
0.474 < q nK {T = 0.7) < 0.637, (30) 
0.368 << ?EA (7 , = 0.8) < 0.556. 

Even with our unprecedentedly long simulations, we are still at the threshold of being able 
to compute ijea- 



5.2 Energy density and scaling exponents 

The time decay of the energy density offers further insights into the connection between 
statics and dynamics (see Sect. [6j. At all temperatures the L = 80 energy data are well 
described by a power-law decay 

e(t w )-e„=AC {T) , (31) 

with e„ the asymptotic (equilibrium) energy value. At the critical temperature T c , general 
scaling arguments relate the decay exponent e to the dimension of the energy operator d e = 
(Dv — l)/v and the critical dynamic exponent z 1 45 ] : 

e(r w ;r c )- eoo;7 - c =Af w £/ " /jc . (32) 

We could then in principle extract e (T c = 1.1) and compare with the best estimate available 
(considering v = 2.45 by Hasenbusch et al. [39] and z c = z{T c ) = 6.86 [33] we expect 
e(r c ) ~ 0.378). Unfortunately we do not have enough statistics at the critical temperature 
to allow for fair fits and statistical accuracy in the determination of the decay exponent, 
but we can avoid such difficulty. Assuming critical dynamics in all the Spin Glass phase 
(T < T c ) such relation could be extended and tested at all simulated temperatures, provided 
a T-dependent dynamic exponent z{T) is considered (see Sect. [3] above). In addition, given 
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T 




A 


e 


0.6 


-1.77862(10) 


0.122(11) 


0.193(6) 


0.7 


-1.77084(7) 


0.148(15) 


0.236(7) 


0.8 


-1.75947(7) 


0.155(16) 


0.268(7) 


0.9 


-1.74429(13) 


0.19(1) 


0.312(5) 



Table 5 Results for the power-law decay parameters, Eq. 43 1 1 - 



an inversely proportional dependence of z(T) on the temperature, we expect that s(T) °c T. 
Once the exponent has been determined for some values of T < T c , we can extrapolate to 
determine e(T~), The exponent may be in principle discontinuous at the critical point, so 
some analysis has been performed also at T = 1.15, slightly higher than T c , to have a hint 
on the possible value of the right limit e(r c + ). 

At all temperatures a fit to power law Eq. Oil ) is quite sensitive to the fitting window. 
This is not surprising as the early dynamics may be very different from the asymptotic be- 
havior. In addition, we are aware that our data may suffer from important finite size effects 
when the coherence length grows up over some fraction of the system size (see [33 1). The 
limit in f w beyond which data for some observable cannot be considered in the thermo- 
dynamic limit depends on temperature and on both the observable and the precision with 
which it can be measured (which for the energy is very high, of order one part in 10 5 ). At 
each temperature, the upper limit in the fitting window should then in principle depend on 
when, in Monte Carlo time, we start experiencing finite size effects in the energy. We can 
check this by comparing with the energy decay of systems of size L = 40. As an example, 
in Fig. [16] we report the difference ez.=4o(f w ) — £l=8o('w) at T = 0.8 as function of f w : no 
appreciable deviations appear in the whole simulated range, allowing us to extend our fit to 
the largest available time. We are then left with adjusting the lower limit f m j n by checking 
that the fit parameters e„, A and e go to stable values. 

Results for e as function of f m j n at all simulated temperatures are depicted in Fig. [T7] 
left. At low temperatures it has been possible to probe for f m ; n ranging in several decades. 
At higher T, it is not possible with our data to have fair fits at larger r m ; n values, as very 




Fig. 16 No appreciable dependence of the decay e(t w ) on the system size, for lattices L = 40 and L = 80 at 
T = 0.8. Inset: the difference between the two as a function of Monte Carlo time. 
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Fig. 17 Left: dependence of e on the fitting window. For each choice of the lower limit t m i u we fit for all 
points t w > t mm . We do not assume any upper limit in the fitting window (see text). At lower temperatures we 
find ; m j n -independent parameters, with clear plateaus showing up. For each curve the plateau value is taken 
considering the midpoint in the last longer series of an even number of points lying on an horizontal line 
within error bars. At T > 1.1 no clear plateau appears, while large fluctuations spoil the fits. The horizontal 
line is a prediction for e(T c ) considering the value of the critical exponent V of 1391 and z c = z(T c ) = 6.86 
from Table|2] Right: the plateau values for e(7") for simulated values of T < 0.9 as function of T (empty 
circles, red color online). A proportionality law e(T) = cT works well, and allows one to extrapolate e up 
to T c = 1.1 (the upper full circle, blue color online). The horizontal line is the best estimate by taking the V 
value of 1 39 1 and Zc from Table|2] 



large fluctuations arise, thus hiding any possible plateau (error bars are from a jackknife 
procedure). 

At T = 0.6, 0.7, 0.8 we considered the parameters converged in the last five points in 
Fig. [TT] left; at T = 0.9 the last three points converged within error bars. In any of these 
intervals, we take the mid-point as the plateau values, and report them in the plot of Fig. [771 
right (corresponding values of all parameters are reported in Table|5]l. The data are well rep- 
resented by a temperature-proportional law, with ^ 2 /d.o.f. = 4.8/3. The extrapolated value 
of £(7^) is 0.373(5) that coincides within errors with the best estimate reported above. This 
is a quite good a posteriori confirmation that the procedure described is robust, leastwise, 
within the precision attainable. From Fig. [17] left, we see that at T = 1.1, 1.15 we were 
not able to identify a clear plateau and estimate e directly, still there is a trend towards the 
expected critical value, supporting that e(T) should be linear up to the critical point at least. 

From the values of e(T) = (3v-l)/(vz(r))0we obtain z(0.6) = 13.4(3), z(0.7) = 
1 1.0(3) and z(0.8) = 9.7(3) in agreement with the results of Tabled 

Alternatively we can link these numerical results with an analytical one obtained by 
Franz et al. 1461 (see also 1481 ). In this reference, the contribution of the interface was 
computed in the framework of Replica Symmetry Breaking. Assuming the power law be- 
havior of the coherence length, one would expect that e(T) = 2.5/z(T). So, inputting our 
values for z(T) we obtain e(T) = 0.33J, which is compatible with our numerical finding 
(e(T) = 0.340(5)7/, see Fig.[[7](left)). 

We conclude this section by testing in nonequilibrium simulations (see also |47]), the 
correction-to-scaling exponent measured by Hasenbusch et al. 1391 CO = 1.0. Our data at 
T c = 1.1 are well described by a double power-law (for f w > 2 6 , CO = 1.0, £ — 0.378) 

e(f w ) = AC (l +ft w 0,/zc ) , (33) 

13 The formula is expected to be valid if the value of the exponents below the critical temperature are 
controlled by their value at T c . 
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giving eoa = 1.70201(2), A = 0.149(2), B = 1.25(5) and x 2 /d-o.f = 104.20/102. 
5.3 The thermoremanent magnetization 

The thermoremanent magnetization of a SG has been known since the 1980s to decay with 
a power law |49 , 50] (deviations to this simple behavior were observed only extremely close 
to T c , namely T > 0.98r c ). As we said in Sect. 12.21 this observable can be identified with 
C(t, f w ) for fixed f w and t ^> r w . Following 1 5 1 j, we have fitted our data to a decay law 

C(r,r„)=A'(f w )+S(f w )f e ('»» . (34) 

The constant term A'(r w ) is justified because for a finite number of samples, the correlation 
function does not go to zero for large times. The very same problem arises in the analysis of 
experimental data 1491 (it was solved by taking a numerical derivative). 

We summarize the results of these fits on Table[6] We fixed f w = 2, 4, 8, 16 and fitted our 
data up to the point where finite size effects appear (t < 2 x 10 10 for T = 0.7, t < 4 x 10 s for 
T = 0.8 and the whole range for T = 0.6). We started our fits at t = 10 6 , but this limit can 
be varied along several orders of magnitude with no change in the fitting parameters or the 
goodness of the fit. As we can see, the exponent c has a slight, but systematic, dependence 
on f w . This tendency was already observed by 1231 . 

Our estimates for c can be compared with experimental values [49 1 



c(T 


= 0.55T C ) s 


t -0.12, 


c(T 


= 0.67r c ) s 


i -0.14, 


c(T 


= 0.78r c ) s 


i -0.17, 



(the error bars on these exponents are small, not much larger than the size of the plotted 
data points in Fig. 3b of [49 1). When compared with the values in Table|6] the experimental 
exponents are similar but slightly higher (let us recall that 0.6 « 0.55r c , 0.7 w 0.64r c and 
0.8 « 0.73r c ). 

A different approach comes from realizing that in experimental work t and r w typically 
differ by at most 4 orders of magnitude, while in our fits they differ by as many as 9 or 10. 



) d(t w ) A'(f w ) x Kg x 1 !^- 

2 TO 5 -0.1525(23) 2.14(5) 2.6(6) 14.7/64 

ig 4 10 6 -0.1495(22) 2.10(5) 2.8(8) 15.5/64 

8 10 6 -0.1459(20) 2.05(4) 2.5(10) 17.4/64 

16 10* -0.1430(19) 2.01(4) 2.4(12) 17.5/64 

2 TO 5 -0.1787(14) 2.067(27) 1.47(25) 23.3/50 

tJ 4 10 6 -0.1765(13) 2.041(26) 1.8(3) 18.4/50 

8 10 6 -0.1733(12) 2.004(25) 1.7(4) 18.9/50 

16 10* -0.1704(12) 1.971(25) 1.6(5) 15.4/50 

2 To 5 -0.210(8) 1.98(9) 1.7(10) 13.9/32 

4 10 6 - 0.212(7) 2.00(8) 2.8(12) 11.1/32 

8 10 6 - 0.208(7) 1.96(8) 3.0(14) 10.8/32 

16 10* -0.205(6) 1.93(7) 3.0(18) 8.43/32 



Table 6 Exponents c and d and parameter A' of equations {34) and 438i . The values of c come from fits, whose 
X 1 we also show. The exponent d comes from d = —cz, where z is the dynamic exponent of equation (22) 
and Table[2] Fits are limited to the time range where the system is free of finite size effects, which accounts 
for the different numbers of degrees of freedom for each temperature. 
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Fig. 18 C(f,f w = 2) as a function of Tlog(t) for three subcritical temperatures. The inset shows that the 
scaling holds only approximately. 



Taking this into account, it is interesting to consider a power law where the fitting window 
is shifted with r w so that 1 < t/t v < 10 (33l . 

C(M W ) = A(/ w ) (1 +f/f w )- 1 / a ('») . (36) 

Using this functional form and extrapolating — l/a(f w ) to a typical experimental time with 
a quadratic fit we obtain 



l/a(*w 


= 100 s) R 


i -0.11, 


T 


= 0.6s 


s0.55r c , 


l/«(fw 


= 100 s) R 


t -0.12, 


T 


= 0.7s 


s0.64r c , 


l/a(f w 


= 100 s) R 


* -0.14, 


T 


= 0.8s 


a0.73r c . 



These extrapolations are slightly above the experimental values of equation P5I >, but if both 
sets of exponents are interpolated with a parabola, the two curves result roughly parallel, 
i.e., our extrapolation error seems temperature-independent. 

Since both the thermoremanent magnetization, Eq. ( 1341 ), and the coherence length are 
well described by a power law, it follows that C(f,f w ) should be a power of | (t + t w ), at least 
for the small values of f w in Table [6] 

C(t,t w )~!;(t + t w )- d . (38) 

Indeed, following the same procedure we used for the exponent a of equation d 1 3 b (see 
Sect. I3.2t we have computed the values of d for f w = 2,4,8, 16 (Tabled, obtaining d « 2. 
Incidentally, let us remark that the coherence length is notoriously difficult to access exper- 
imentally [5 |, while the thermoremanent magnetization is a pretty standard measurement. 
Eq. d38b then appears as an interesting, albeit indirect, way of estimating experimentally an 
effective coherence length. 

Notice that the exponents c(f w ) of Table [6] are roughly, but not exactly, linear in T. 
This suggests that the thermoremanent magnetization should be a temperature-independent 
function of T\og{t) 1501 . Fig.[T8]shows that this is only an approximate claim. 



2S 



T 


tw 


i n 


e(t w ) 


/Ow) 


^/d.o.f. 




2 


10 J 


-0.236(7) 


0.873(9) 


52.2/104 






10'' 


-0.30(6) 


0.82(5) 


13.9/64 




4 


ur 


-0.203(6) 


0.909(8) 


47.9/104 


u.o 




10 6 


-0.25(4) 


0.85(4) 


13.5/64 


8 


10 3 


-0.176(4) 


0.943(7) 


41.5/104 






10 6 


-0.21(3) 


0.90(4) 


13.9/64 




i /; 

lo 


ur 


-0.158(4) 


0.968(7) 


38.1/104 






10 h 


-0.19(3) 


0.92(4) 


15.3/64 




2 


10 J 


-0.263(4) 


0.890(4) 


43.0/90 






10 6 


-0.32(3) 


0.84(3) 


14.4/50 




4 


10 3 


-0.230(3) 


0.921(4) 


71.9/90 


7 




10 b 


-0.29(3) 


0.862(25) 


12.8/50 


8 


ur 


-0.2003(23) 


0.955(3) 


94.1/90 






10 6 


-0.253(23) 


0.895(23) 


13.0/50 




10 


10 3 


-0.1768(20) 


0.985(3) 


138/90 






10 6 


-0.226(19) 


0.921(22) 


10.6/50 




2 


10 J 


-0.302(16) 


0.891(16) 


45.1/72 






10'' 


-0.5(3) 


0.77(15) 


14.3/32 




4 


ur 


-0.257(12) 


0.934(14) 


63.6/72 


0.8 




10 6 


-0.6(4) 


0.71(17) 


11.4/32 


8 


10 3 


-0.223(10) 


0.970(13) 


69.8/72 






10 6 


-0.49(24) 


0.76(12) 


11.1/32 




16 


ur 


-0.192(8) 


1.008(12) 


65.9/72 






10 s 


-0.40(19) 


0.81(12) 


8.49/32 



Table 7 Parameters of a fit to Eq. {39), offering an alternative description of the thermoremanent magnetiza- 
tion. For each temperature, the maximum time included in the fit was such that ^ (f max ) = 10. 



The incompatibility of the values in Eqs. 05) and (37) . together with the fact that we 
needed the constant term A'(f w ) in Q4b . suggests that there probably exists a systematic 
error in the power law fits of Table ©Lij One may consider an alternative description, 

C(f,f w )=A"(f w )exp[ e (f w )(logf) /( ' w) ] , (39) 

that would reproduce a power law if /(r w ) = 1. Eq. 439b should not be confused with the 
stretched exponential, discarded in reference |49] for all T but the closest to T c . We show 
our results for this fit in Table [JJ From the point of view of a test, the two descriptions, 
( 1341 ) and (39) . are equally good. The fit parameters are remarkably stable with variations of 
the fitting window. As we see, the values of f(t w ) are very close to, but incompatible with, 
1 (at least for the lowest temperatures). 

Let us conclude this section by considering again the effects of statistical data correlation 
on fit parameters. Specifically, let us consider the exponent cc(t w ) of equation (36) , see 
Fig.[H]and reference [33]. The alert reader will notice strange wiggles, large as compared 
with the error bars, which are specially prominent for the fit with 63 samples. The reason 
is that data for different t and f w in this fit are even more correlated than usual. In fact, we 
have sampled the function C(f,f w ) for t and ? w integer approximations to 2 ,//4 , i = 0, 1,... 
For each r w , we perform the fit and extract a(f w ) for r w < t < 10r w . In other words, for each 
t w we used 14 values of t. Recall that the spin configurations involved are the one at time 
t w , and the 14 spin configurations at succeeding times t + t w . Given our choice for integer t 

14 Note, however, that the smallest value of C(t,t w = 2) that we reach is ~ 0.013, pretty large as compared 
with A'(f w ). 
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Fig. 19 Exponent <x(t„) defined in equation j36\ as a function of t w at T = 0.7, computed for the 63 samples 
of [ 33 ] and for our 768 samples. The right panel shows the values of diagonal y} for both fits. 

and r w , it follows that, for the next r w , in the computation of C we will use 13 out of the 14 
spin configurations used at the earliest f w . For the second-next f w , the number of repeated 
configurations will be 12, and so forth. This is the origin of the dramatic data correlation: 
the very same spin configurations are being used for consecutive times. Notice that when the 
statistics is increased to 768 samples, the period of these oscillations does not change, but 
their amplitude decreases. Also in Fig.|T9j right panel, we show the value of % 2 ld.o.f. for the 
fits with 63 and 768 samples. In both cases, % 2 decreases strongly with f w , but, while the fit 
with 63 samples was perfectly reasonable from r w ~ 10 4 , that with 768 samples is not good 
until t w ~ 10 s . The increased accuracy reveals systematic deviations from equation ( 136b . 
Nevertheless, the estimate of a(f w ) for both fits is compatible in a much wider range. 



5-4 Cu n k as a function of C 2 

Studying Qink as a function of C 2 for fixed f w we can monitor the scaling of the active 
domains' surface area. This scaling has an exact correlate in equilibrium studies [52}, where 
one considers the probability distribution function of 

Gnnk = ^ £ aj V 2 W% (2) • (40) 

The analogue of C\i,± as a function of C 2 is the conditional expectation value (Qiink|<7 2 )- 
Indeed, in 1331 . we quantitatively compared our dynamical results with the equilibrium 
( Clink | <7 2 ) of 1521 . It was found that a convenient time-length dictionary could be estab- 
lished in such a way that the equilibrium results for finite lattices reproduced the nonequilib- 
rium results for finite £,(t w ). For instance, for T = 0.7, when L/| (r w ) w 3.7, Cn^C 2 ,f w ) ~ 
(Giink|<7 2 )L- 

In 1331 it was found that Cii n k is a smooth increasing function of C 2 at least up to ex- 
perimental scales. On the other hand, for systems undergoing coarsening dynamics (e.g., a 
disordered ferromagnet), Cy n ± tends to a constant C 2 -independent value whenever C 2 < #| A . 
Let us briefly justify these expectations. 

On the one hand, since in the RSB scenario the coherent domains are not compact ob- 
jects, one would expect C]j n k to have the same aging properties as C 2 , that is, dCi; n k/dC 2 
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should not vanish. This is the nonequilibrium analogue of the overlap equivalence prop- 
erty 1 52 1 . For instance, in the Sherrington-Kirkpatrick model it is straightforward to show 
that dink = C 2 . 

On the other hand, to find the scaling for a coarsening image of compact active droplets 
we need a more elaborate argument. We consider a large droplet of size ^ (f +f w ) at time 
t +r w that, at time r w , was made of Nq smaller droplets of size <§(f w )- The number of spins 
in the boundary of a droplet at time r w scales as <§(f w )° s . Typically, D s = D — 1, but one 
may have D - 1 < D s < D (6lfl6l (for the TNT model of SG, D - D s w 0.45 for D = 3, 
see 1 55 1 and Palassini and Young in [ 18 1). Of course, Nc scales as Nc ~ [%(t + f w )/<§ (f w )] D - 
The overlap of each of the Nc droplets at time f w with the configuration at time t + f w is 
randomly ±<?ea- Hence, the scaling of C(?,? w ) is, for the region C < q'eaIH 



^(f + fw)/ \%{t + t, 

Now, for the link overlap we expect (C[' nk is the equilibrium expectation value of Qiink) 

C, mk (f,/w) = < k + N c ZDU + l) ■ ^ 

In fact, the decay of Cv m k(f,fw) comes mainly from the contribution of droplets' surface at 
time r w . In particular, for t — > °°, the excess of Ci; n k over Cj? nk is just the probability that a link 
belongs to the surface of a droplet at time f w . Now, considering equation ( 141b we conclude 
that the number of droplets Afc scales with C(t,t w ) as 

Nc- 8 -^-, (43) 

where the function g{C) is continuous, but not necessarily differentiable at C = 0, and where 
g(0) > 0. Combining Eq. (EUl and Eq. (gT] we get 

Ci mk (Mw) = (3L + cUg(C)t; D *- D (t w ) . (44) 

In the above expression C/ ink is a constant. Notice that, in particular, equation d44b implies 
that the derivative of Ci; n k with respect to C 2 goes to zero as f w — > °°. 

We can easily check Eq. ( 144b for the two-dimensional Ising model, where the Onsager 
and Yang solutions provide exact values for Cj] nk and (in this case D — D s = 1). As for 

the growth of the coherence length, it scales as ?4 (see, for instance 11531 ). As expected, see 
Fig.HQl Cii n k tends to a constant function Ciink(C,f w — > °°) = C° nk for C < qEA- As we show 
in Fig. [20] right, the approach to this constant is well described by Eq. J44b . 

As for the Edwards-Anderson model, Ci; n k is seen to be a very smooth function of C 2 1331 , 
so we can easily compute the curve Cii„k(C 2 = 0, f w ) with a linear extrapolation. We can also 
study Ci(r = l,f w ), which is the nonequilibrium disorder average of Qii n k(fw)> Eq. ( 140b . The 
two curves are plotted against | _1 (t w ) and ^ _0 45 (f w ) in Fig.|2TJ In accordance with the pre- 
vious discussion, and in particular with the relation C nn k(C 2 = 0,f w (£)) = (Qimk\q 2 = 0)l, 
both have the same extrapolation to infinite time (they actually collapse on the same curve 
for large times). Of the two, the <§ _1 scaling is more convincing. 



15 Even if Eq. 141 1 is intuitively evident, it can be backed by an explicit computation. From Eq. (6-11) 
in [ 54 1, one easily shows for the Ising ferromagnet that the spin-spin correlation function takes the form of a 
series C(t , t w ) = a iy + a 2 y 2 + . . . , where y = % °/ 2 (( + f w ) % °/ 2 (f w ) / (f + t w ) + ^ (f w )] D / 2 . Eq. gD follows 
when^(? + r w )>^(f w ). 
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c c 

Fig. 20 Left: Qi n k as a function of C for the two-dimensional ferromagnetic Ising model at T = 0.66T C . 
Results obtained for an L = 4096 lattice (the results were averaged over 20 trajectories). Right: Numerical 
check of Eq. (44) for the data on the left panel. The vertical line is at C = <?ea- We see that, for large f w and 
C<q E A, CjnkfCfa) - scales as i; -1 . 
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Fig. 21 Extrapolation of Cii n k (C ,f w ) to C = 0, together with Ct(r = 0,f w ), against 7 (( w ) (Left) and 
^f2' 45 ('w) (Right) (plot for our 96 samples at T = 0.6). 



Equation | |44| > suggests plotting dCii n i c (C 2 ,f w )/dC 2 (see (33ll for details) against <§ _1 (f w ) 
(Fig. [22] left) and against ^ _0,45 (f w ) (Fig. l22l right). It is important to choose a value C 2 
of C 2 smaller than g| A but not too small, because otherwise the numerical estimate for 
the derivative would be unreliable. We have used the most pessimistically small estimates 
of q\ A in Table [4] The two representations are linear within our errors. However, while a 
£,~ (f w ) scaling compatible with standard coarsening seems falsified (the extrapolation to 
infinite time is well above zero), a ii;~ - 45 (? w ) scaling towards zero is compatible with our 
data. 

From Fig. [21] and Eq. J441) . we conclude that the difference C\\ n ^{C 2 ,f w ) — C^{r = l,f w ) 
should vanish as <§ Ds_Z5 (f w ) in a coarsening system. We show this difference in Fig. [23] for 
the same fixed value C 2 we used for the derivative. As we see, the extrapolation as <§ _1 is 
smooth and positive. On the other hand, the ~° 45 scaling could only be possible if our 
whole simulation were in a pre-asymptotic regime. 

This discussion notwithstanding, two comments are in order. First, equation ( 144b relies 
on equation d41 b . which is disproved by the values of exponent d in Table [6] (we obtain 
d w 2, rather than d = 3/2). Second, we mark by crosses in Fig. [22] the experimentally 
relevant scale: the derivative is certainly nonvanishing there in either case. 
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Fig. 22 Derivative of Cii„k(C 2 ,f„) with respect to C 2 versus £~'(f w ) (Left) and versus | _0 ' 45 (f w ) (Right) 
for three subcritical temperatures (T = 0.6,0.7,0.8, from top to bottom). Lines are linear least squares fits. 
We mark by crosses our extrapolations for the experimental scale of £, (f w = 100 s). The curves are plotted for 
a fixed value C 2 = Cj, chosen to be just below our lower bound for ^ea at each temperature from Eq. )30> . 




Fig. 23 Difference Ci ink (C 2 ,f„) - C 4 (r = l,f w ) for the same C\ of Fig.|22] The curves are plotted against 
£ _1 ( ( w) (Left) and against i;~ ' 45 (f w ) (Right). An extrapolation to zero seems unlikely even for the £~ a45 
case. 



6 Scaling of the dynamical coherence length 

As we have shown in Table [2] our data for the coherence length is very well represented 
by a power law in f w . We found an exponent z{T) roughly linear in T^ 1 . Some theoretical 
grounds for this behavior can be found in I51II56II . It also appears in the numerical simulation 
of the Sherrington-Kirkpatrick model 1 57 ]. Nevertheless, an alternative interpretation has 
been suggested 1 12 1. We now reanalyze our data under this light. 



6.1 Mixed scaling 

The Saclay group proposed 1121 . see also (9), a mixed scaling for the dynamical coherence 
length, which assumes both critical behavior and activated dynamics in a wide range of 
temperatures in the glassy region: 

f w ~T ^exp(^^V (45) 
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Fig. 24 Left: Function G{t w ,T) denned in Eq. {46) versus ( w , for T = T c = 1.1, T = 0.8, T = 0.7 and 
T = 0.6. Our estimates for the plateau (see text) are indicated with horizontal lines. The continuous curves 
are G(t w ,T) as computed from {48), fixing Tq(T) = To(T c ), using z(T) computed in the power law fits of 
Table|2] At T = 0.8 we also show our data for L = 40. Right: As in left panel, but now we allow To to depend 
on T, Eq. {49) . 



where To is the microscopical time associated to the dynamics; z c is the dynamical critical 
exponent computed at the critical point; \y is the exponent that takes the free energy barriers 
into account (from the dynamical point of view) and Y(T) = Fq(1 — T/T C )V V , with the v 
exponent being the static critical exponent linked to the coherence length. Near the critical 
point Y(T) — > and the power law critical dynamics is recovered. 

To asses the validity of the mixed scaling hypothesis, we consider the following func- 
tion fl2l : 

'log(r w /T ) -z c log<§(f w ,r) 



G(/ W ,r) 



1 

y/v 



S*T e /T 

Equation ( 145b would imply that G(/ w , T) is a f w -independent function of temperature, 



G(^w,7 , ) = Go 



T 

1 



(46) 



(47) 



where Go = {Yq / T^ 1 1 Both the Ising and Heisenberg experimental sample behave 
consistently with this expectation. In this study, the parameters were taken to be z c = 5 and 
V = 1.5 (those of AgMn, a Heisenberg SG) fT2l . 

Notice that G(r w , T) would be exactly zero if followed a pure power law with z and To 
fixed to their critical temperature values, 



1 Vz(Tc) 



(48) 



In the case of a power law with parameters T and z different from those computed at the 
critical point, we expect that the function G(t w ,T) should be zero only for t w — ► °° . In order 
to avoid a painful and somewhat arbitrary fit, we take the relevant parameters in Eq. ( |46t 
from the literature (z and T are taken at T c ): T c = 1.109 (39), y ~ 0.7 (23), z c = 6.86 (33) 
and v = 2.45 (39). 

In Fig.[24]left, we show the function G(? w , T) for four values of the temperature (includ- 
ing T c ). Although G(f w , T) is not ^-independent, it plateaus at a value G (T) for long times. 



16 Ag:Mn at 2.5% (Heisenberg like), CdCr1.7IN0.3S4 (also Heisenberg like) and Feo.sMno.sTiOj (Ising 
like). 
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Fig. 26 Coherence length 2(*w) as a function of riog( w , for three subcritical temperatures. Even if the 
three curves are not equal within errors, the overall scaling is suggestive. 



This is especially clear for T = 0.6 and T — 0.8, while at T c we expect it to be compatible 
with zero. For T = 0.7 the plateau is not well defined, but we estimate it as the average of 
the last points^- As we show in Fig. [25] G a (T) behaves consistently with Eq. 047 | l. Hence, 
we are in a time region where experimental results for G(/ w , T) are reproduced. 
However, we will remark the following points: 

1. The departure of the curve for T = 0.8 from the plateau in Fig.[24]is a finite size effect, 
as shown explicitly by our L = 40 data (see also 1331 ). A similar, though milder, effect 
afflicts the data for T = 0.7 for r w > 2.2 x 10 10 . The same effect is very clear at T c , for 
even shorter times. 



In order to assign error bars to the plateau values, we consider only the biggest contributions, those from 
the uncertainties in z c and in To(r c ). 
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2. Were (f w ) a power law, the curve for G(t w ,T) would be decreasing (for larger times). 
On the other hand, we know that finite size effects cause it to increase. The plateau may 
well be an artificial combination of these two effects. 

To clear up this problem, we also show in Fig.[24j right, the value of G(/ w , T) computed 
using the coherence length and Tq(T) from a fit to 



1 l/z(T) 

(49) 



Note that, at variance with Eq.d48b. we now allow To to depend on T. These fits were reported 
in Table H Specifically, we obtained T Q (T = 0.6) = 0.008(3), T (T = 0.7) = 0.030(8), 
T (r = 0.8) = 0.17(4) and T (r c ) = 0.58(13). 

Although the physical meanings of Eqs. l !46l > and l !49l > are quite different, the two of 
them account fairly well for our data. 

However, Eq. fl49b describes well the data in nearly the whole range of times (actually 
for <^ > 3) while Eq. d46l > describes the data only in the region where the right hand side 
of Eq. 046 l l, computed with Eq. l ]48t , is nearly constant. Moreover the data for £, (r w ) nearly 
collapse if we use the variable Tlogr w , see Fig. [26] This collapse is not compatible with 
Eq. i45\ . Henceforth, the activated scaling hypothesis, Eq. i45l , implies that our data are 
entirely in a pre-asymptotic regime. We note nevertheless that extrapolating our data with 
Eq. I l49t to the relevant experimetal scale (r w = 10 14 or 100 seconds) produces fairly sensible 
results [33 1. 



7 Conclusions 

Using the dedicated computer Janus, we have studied the nonequilibrium dynamics of the 
Ising spin glass for times spanning eleven orders of magnitude. We have looked into quanti- 
ties not considered in our previous work 1 33 1 and extended the simulations described therein 
by considering more temperatures and vastly enlarging the number of samples for T = 0.7. 
The emerging picture is that of non-coarsening dynamics. 

We have performed an extensive investigation of heterogeneous dynamics, by consid- 
ering the two-time, two-site correlation function C2+i{r ,f ,r w ) . We have obtained the first 
reliable determination of the nonequilibrium correlation length and the exponent for the al- 
gebraic decay of Ci+2- When t is much smaller than f w , the correlation length reaches a 
r w -independent limit. On the other hand, for t much larger than r w , the correlation length 
scales as the coherence length. Thus, it might be sensible to exchange the role of both length 
scales in the study of structural glasses 1421 . where a notion of a coherence length is lacking. 

Crucial to the above findings has been our integral determinations of characteristic 
length scales. We have also used them to obtain the coherence length and to study the 
replicon mode. Indeed, the exponent a in Eq. dl3l l is definitively nonvanishing, and hence 
incompatible with the droplets picture. We have also considered nonequilibrium overlap 
equivalence, with the help of the coherence length. 

We have used both the coherence and correlation lengths to obtain safe bounds for the 
Edwards-Anderson order parameter below the critical temperature. 

18 We have found a monotonic (decreasing) behavior for the microscopic time just as in Ising samples 
(Feo.5Mno 5T1O3) [58]. However, Heisenberg samples have no such clear pattern: CdCi'i 7IN0.3S4 shows 
decreasing monotonic behavior but Ag:Mn at 2.5% and Cu:Mn at 6% [5 ] present an increasing monotonic 
one. 
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As for the thermoremanent magnetization, good agreement with experimental determi- 
nations of the temperature-dependent decay exponents have been obtained. A potentially 
useful observation for experimental work is that the thermoremanent magnetization scales 
with the coherence length, which is much harder to measure. We also observed that a non 
power law function could fit the thermoremanent magnetization equally well. 

The energy relaxation is well described by a power law. The exponents displayed a 
nearly linear dependence on temperature. It has been possible to extrapolate to the critical 
point, finding results in nice agreement with the latest determinations I33II39I . 

We have shown that the link overlap correlation function Cii n k offers a wealth of infor- 
mation on interphase behavior. Our results have been equally compatible with the droplets 
and RSB pictures. However, in the droplets picture the scaling with the coherence length of 
the thermoremanent magnetization is incompatible with our data. Furthermore, irrespective 
of what happens for infinite time, the variation of Ci; n k with C 2 is nontrivial at experimen- 
tal time scales. This means that the physical view conveyed by the RSB theory is a better 
representation of the physics at the scale of a few hours. 

We have critically examined the time growth of the coherence length, comparing critical 
and activated dynamics. We have found that both theories describe its behavior equally well. 

Finally, we have taken the occasion to give full details of our analysis methods, some of 
which are quite new. 
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